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Chains of magnetic atoms placed on the surface of an s-wave superconductor with large spin-orbit coupling 
provide a promising platform for the realization of topological superconducting states characterized by the 
presence of Majorana zero-energy modes. In this work we study the properties of one-dimensional chains of 
Yu-Shiba-Rusinov states induced by magnetic impurities using a realistic model for the magnetic atoms that 
includes the presence of multiple scattering channels. These channels are mixed by spin-orbit coupling and, via 
the hybridization of the Yu-Shiba-Rusinov states at different sites of the chain, result in a multi-band structure 
for the chain. We obtain the topological phase diagram for such band structure and show that the inclusion of 
higher bands can greatly enlarge the phase space for the realization of topological states. 

PACS numbers: 73.20.Hb, 74.78.-w, 75.70.Tj, 


I. INTRODUCTION 

The search for Majorana zero modes (Majoranas) in con¬ 
densed matter systems has been an active and exciting pur¬ 
suit [1-4]. The reason for much excitement is due to the 
theoretical prediction that these modes manifest non-Abelian 
quantum statistics [5-8], and, as such, would open the possi¬ 
bility to realize topological quantum computing [9-12]. Cur¬ 
rently, errors caused by the decoherence of the quantum states 
used to encode the data constitute the biggest fundamental 
obstacle for the realization of a scalable quantum computer. 
In topological quantum computing this obstacle is overcome 
by the topological protection of the quantum states used to 
encode the information. Most platforms for realizing topo¬ 
logical phases of matter supporting Ising anyons (i.e. exotic 
defects binding Majorana zero modes) involve superconduct¬ 
ing heterostructures [13-33]. Recently, several works [34- 
51 ] have proposed that a chain of magnetic atoms placed on 
the surface of a superconductor can be in a robust topological 
phase characterized by the presence of Majorana modes lo¬ 
cated at its ends. In addition, recent experimental results [52] 
have shown that a chain of Fe atoms placed on the surface 
of superconducting Pb exhibits a zero-bias peak localized at 
its ends, consistent with the presence of Majorana modes. A 
very recent preprint also presents experimental results for a 
similar system [53]. A full understanding of the experiment 
presented in Refs. 52 and 53 is still being developed, for exam¬ 
ple, the height of the zero-bias peak is a small fraction of the 
predicted universal value of 2e^/h [54, 55], a fact that could 
be attributed, for example, to finite temperature broadening or 
disorder. 

The potential of a chain of magnetic impurities, placed on 
the surface of a superconductor, to be in a robust topological 
phase, calls for a thorough theoretical understanding of this 
system. So far most of the works have assumed one bound 
state per magnetic atom, corresponding to the zeroth angu¬ 
lar momentum channel (/ = 0) of Yu-Shiba-Rusinov(YSR) 
states [56-58]. However, realistic adatoms are expected to in¬ 
duce several bound states corresponding to different angular 
momentum scattering channels (i.e. I = 0,zbl, etc.). It has 


been shown experimentally that partial waves beyond s-wave 
are essential to explain the energy spectrum of magnetic atoms 
such as Mn, Cr and Fe [59-61]. Furthermore, the YSR states 
originating from the / = 0 channel (s-wave) are not always the 
lowest energy eigenstates. It is therefore necessary to under¬ 
stand the interplay of different angular momentum channels, 
in particular in the presence of significant Rashba spin-orbit 
coupling(SOC) [62]. The Rashba SOC term in the effective 
low-energy Hamiltonian describing the fermionic degrees of 
freedom is expected to be present at the surface due to broken 
inversion symmetry [63]. Moreover, the presence of SOC is 
required in order to have a stable topological phase for a fer- 
romagnetically ordered chain [42, 52] placed on the surface of 
an s-wave superconductor. 

In this work we study a realistic model for a chain of mag¬ 
netic impurities placed on the surface of an s-wave supercon¬ 
ductor, see Fig. 1, that takes into account multiple scattering 
channels for the adatoms and the presence of SOC. Treating 
the magnetic impurities classically we study the properties 
of the bands formed by the hybridization of the YSR states 
bound to the different adatoms forming the chain. The multi¬ 
channel treatment of the scattering potential of a single impu¬ 
rity implies that, for the chain, we obtain a multiband model. 
In the remainder, to simplify the presentation and to be able 
to clearly point out the main qualitative features of the multi¬ 
band structure resulting from the hybridization of multichan¬ 
nel YSR states at different impurity sites, we assume that the 
YSR states corresponding to different values of |/| are well 
separated in energy. In this limit we immediately notice a fun¬ 
damental difference between the band, s-band, formed from 
the |/| = 0 states and the bands formed from |/| > 0 states: 
the s-band can be assumed to be well separated from all the 
other bands and therefore can be treated effectively as a single 
isolated band; however, given the degeneracy, in the limit of 
no SOC, of the +|/|, —1/| states, we have that the two bands 
formed from +1 /1, and — | /1 states are always very close in en¬ 
ergy and therefore that for bands formed from |/| > 0 states, 
a multiband treatment is necessary. Henceforth, we limit our¬ 
selves to the case in which max(|/|) = 1, and consider the 
minimal model that captures the aforementioned features. 
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The topological properties of the chain are determined by 
the band (bands) that is (are) closest to the midgap energy of 
the superconductor. To exemplify the main properties of the 
chain of multichannel magnetic impurities we consider two 
limits: (i) the limit in which the I = 0 YSR states are closest to 
the midgap energy of the superconductor so that for the chain 
the resulting s-band is also the closest to the midgap region, 
see Fig. 1 (b); we call this case “deep s-band” limit; (ii) the 
limit in which the |/| = 1 YSR states are closest to the midgap 
energy of the superconductor so that for the chain the resulting 
bands, “p-bands”, are also the closest to the midgap region see 
Fig. 1 (c); we call this case “deep p-band” limit. We obtain 
the topological phase diagram for both the deep s-band and 
the deep p-band limit. The presence in the deep p-band limit 
of two bands close in energy would suggest that in this limit 
the topological phase could be strongly suppressed. Contrary 
to this naive expectation we find that in the deep p-band limit 
the phase space in which the topological phase is present can 
be even larger than in the deep 5-band limit. We also find 
that in the multichannel case the presence of SOC leads to the 
dependence of the chemical potential on the direction of the 
chain magnetization, which, in principle, allows one to tune 
between topological and non-topological phases. This is an 
important feature for braiding Majoranas [64]. 

The paper is organized as follows. In Sec. II we introduce 
our model and explain the general framework for the calcula¬ 
tion. In Secs. Ill and IV, we derive an effective Hamiltonian 
and calculate the topological phase diagram as well as quasi¬ 
particle gap for the deep s- and p-band limits, respectively. 
In Sec.V we discuss the qualitative difference between those 
two limits and compare their topological phase diagrams. The 
technical details are presented in the Appendices. 



(b) Deep s-band limit (c) Deep p-band limit 




-A 






FIG. 1. (Color online) Schematic setup of the proposed structure 
supporting multichannel Yu-Shiba-Rusinov chain. Atoms with high 
magnetic moments form a ferromagnetic chain on the surface of thin 
film SC [42, 52]. Each magnetic atom creates multiple Yu-Shiba- 
Rusinov states with distinct angular momentum quantum number 1. 
Here we consider I = —1, 0, and 1 states which forms three bands 
in a quasi one-dimensional system. In the topological phase, there is 
an odd number of Majorana zero-energy modes at the opposite ends 
of the chain. 


II. THEORETICAL MODEL FOR MULTICHANNEL YSR 
CHAIN 

We consider a chain of magnetic impurities separated by 
a distance a and placed on top of an 5-wave superconductor 
with Rashba spin-orbit coupling (SOC). The corresponding 
Hamiltonian describing an effectively two-dimensional super¬ 
conducting film with Rashba SOC reads (h = 1): 

Hsc = [^kCTo + Of (o- X k) • z] (1) 

where k = ky) the electron momentum, = 

k^!2m — iJ. with /i the chemical potential and m the effec¬ 
tive electron mass, a is the strength of the Rashba SOC, z is 
the unit vector normal to the plane, see Fig.l, and A is the 
superconducting gap. The Pauli matrices and operate in 
spin and particle-hole space, respectively. The presence of the 
magnetic impurities is taken into account via the Hamiltonian 

^fimp = ^ y,- (r - R, ) = - y] J(r - R,.) (S, • <7) tq ( 2 ) 

3 3 

where Vj are the individual impurity magnetic potentials with 
Rj, Sj and J being the impurity position, its classical spin 


and its exchange coupling to the host quasiparticles, respec¬ 
tively. In order to find the band structure of the chain of YSR 
states induced by the magnetic impurities, we need to solve 
the Schrodinger equation: 


[^fsc + ^fimp(r)]t/>(r)=i;t/>(r). (3) 


Here ^/^(r) is the Nambu spinor 

^ ) • The spectrum of the 

subgap states is determined by the pole of the T-matrix [65]. 
In terms of the Green’s function for the superconductor 
G = [E — the eigenvalue problem for the subgap 

states is given by (1 — GHimp)'ip{T^) = 0. We note that the 
Green’s function of an s-wave superconductor in the presence 
of SOC has both even- and odd-parity components [66], 
and can be written as G{k;E) = ^iGi{k; where 

k = \k\ and = aictdinkx/ky. In the presence of Rashba 
SOC, the Gi are non-zero only for / = —1, 0,1 and are given 
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by: 


G-i{k-,E) 


Go{k-,E) 


Gi{t,E) 


1 . •>X Ea+To + Aa+Tx + ex{k)a+Tz 

2^^^ ^ ’ E^-el{k)-A^ 

(4) 

1 EaoTo + AaoTx + s\{k)aoTz 

2^^ E^-sl{k)-A^ 

1 Ect^tq + Aa_T^ + ex{k)a-T^ 

2^^* E^-sl{k)-A^ 

( 6 ) 


where a± = {a^ ± iuy) /2, and e\{k) = fc^/2m — \ak — ji 
are the dispersions of the helical bands with A = ±. 

Transforming to momentum space, Eq. (3) can then be 
rewritten as: 


^i(k) = 

3 


(7) 


where 'ipi(k) is the Fourier transform of the wave function 
centered at . Assuming that at the Fermi surface the scat¬ 
tering potential is weakly dependent on p = |p| we have 
V^ (k — p) ~ Vji^k — ^p), where 




-JiOk - Op)Sj ■ a 
0 


-J(0p 


0 ) 

- ^k)S,■a J 


( 8 ) 

The magnetic potential can be decomposed into inde¬ 
pendent angular momentum channels: Vj{0\^ ~ ^p) = 

so that 


V.' 


hi 


Ji^j ' ^ 
0 


0 

—J_/Sj 



(9) 


with Ji the angular momentum components of J(6>k — ^p). 
Since Vj{0\^ ~ ^p) is Hermitian and an even function of 
— ^p, we have Ji = J-i and Vj^i = —Ji {Sj ■ a) tq. We as¬ 
sume that I = 0, the s channel, and I = ±1, the p channels, are 
the dominant scattering channels: Vj = X]z=o ±i ^ 3 ,i- For the 
impurity chain arranged along the x direction in a ferromag¬ 
netic pattern, i.e., Sj = (cos(p) sin^, sin(p) sin^, cos(^)), 
(illustrated in Fig. 1(a)), Eq. (7) in the angular mode repre¬ 
sentation takes the form 


hi' 


where an overline indicates an integration over momentum 
amplitude: = f with ipi^i the angular mo¬ 

mentum components of '0i(k) ('0i(k) = 

2in& = j ^dk^^ ^ikxij cosei^^-i{i-i')ei^Qfj^.^ 

Xij = Xi — Xj. Using Eq. (4)-(6) one can rewrite in 

terms of the dimensionless integrals: 

gi/CA i£)x COS 0k g7n0k ^ 

^2 - ^2 - A2 ’ 

( 11 ) 




Kn,\{x\E) = 


Nx 


27r2Ar, 


,j:m 


^ ^ik\{£)x COS 9\^^in0\f^^ 

de- 

-D 


p;2 _ g2 _ 


( 12 ) 

where D is an energy cut-off, kx{£) = kp^x E £/vf,x with 
kp^x — kp (Vl “h + Aof), vfx — 'zzfa/I E kp = 
vf = kpIm. a = am/kF is the dimensionless SOC 


coupling. Nx = 


1 + A 


\/T+^ 


is the density of states of 


the A helical band at the Fermi level in the normal state, and 
Nf = {Np + N-)/2. The analytic results for the above in¬ 
tegrals in the limit D ^ oo and for Qi^_i,{E) are presented 
in Appendix A and B. For convenience, we define the dimen¬ 
sionless exchange couplings Ji = Ji\S\7tNf which will be 
used henceforth. 

It is convenient to rewrite Eq. (10) in the following form: 


^M^^{E)^j = 0 


(13) 


where = {'ip^ is a. 12 dimensional spinor, 

and the matrix M.E[E) is defined as = SijSi^i^ — 

Qy_i,{E)Vi'. Here the local part of the matrix de¬ 

termines the YSR spectrum of a single magnetic atom [62] 
whereas the non-local part describes the hybridization 
between YSR states induced by the magnetic atoms at i and j 
sites. For an equally spaced magnetic atom chain with dis¬ 
tance a between the two nearest atoms, this hybridization 
leads to the formation of the YSR bands. In the limit of 
kpa ^ 1, which we consider henceforth, the hopping energy 
scale is proportional to 1/ ^/k^ and, thus, the bandwidth W 
is small, i.e. W ^ A. In this limit, the bands maintain the 
character of the single impurity YSR states and, thus, we re¬ 
fer to them as s or p-bands. Strictly speaking, SOC mixes 
different angular momentum states but, since we assume that 
a <^1, this terminology is justified. 

When s and p bands are well-separated by a gap that is 
much larger than the temperature, see Fig. 1(b) and (c), the 
problem can be considerably simplified by integrating out the 
higher-energy bands (i.e. by taking into account virtual scat¬ 
tering processes to the bands higher in energy). In the follow¬ 
ing, we consider two limiting cases corresponding to the deep 
s- and p-band limits and discuss the corresponding topologi¬ 
cal phase diagrams. We show that these two cases are qualita¬ 
tively different since deep p band limit consists of two bands 
originating from the / = ±1 YSR states. 


III. DEEP S BAND LIMIT 
A. Derivation of the Effective Hamiltonian 

We first consider the deep 5-band limit such that the energy 
of the / = 0 state is close to the midgap, i.e., Jo ^ 1 with the 
on-site energy cq ~ A(1 — Jo) -f 0{a‘^) 0. Provided Ji <C 

Jo the / = 0 band is well separated from the p bands, i.e., the 
bandwidth W C A | Jo — Ji |. After integrating out the I = 
±1 states, we obtain a tight-binding description for the single 
5-band with the p channels taken into account perturbatively. 
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by allowing for virtual transitions through the p channels. This 
can be done by first solving for using Eq. (13) 

i’i,i = + E (14) 

and substituting above expressions into the equation for / = 0 
component and keeping terms up to the linear order in inter¬ 
site coupling, we obtain 

Em:^'W,.o = 0, (15) 

3 

where the matrix M*-^ (£’) is given in the Appendix C. In or¬ 
der to solve Eq. (15) analytically, we expand the local on-site 
matrix to the linear order in E around = 0, assuming that 
^0 ^ 0? 

• E, 

and set E = 0 in the inter-site matrix: 

lim Mf^(E) = 

E^O ^ ^ ^ s V / 

In doing so we ignore terms 0{l/kFCL) <C 1 and 
O . With these approximations, Eq.(15) can be 

written as 

= (16) 

3 

where the local and non-local contributions are 
given by W/ = and Wj = 

M*^-^(0), respectively. The tight-binding 

Hamiltonian 1-Ls{hj) is obtained by projecting Eq. (16) 
onto the local YSR states: ((/:?+, where ip± are 

the eigen spinors of the single-impurity bound states with 
energy Ecq. The local basis can be found by solving the 
single-site equation M^*(E')(/9± = 0 as a special case of Eq. 
(15), where the bound state energies are determined from 
Det [M**(E)] = 0, see Appendix D. 

We first assume that the magnetic atom moments are 
aligned ferromagnetically along a) z- (out-of-plane), b) x- 
(along the chain direction), and c) y- (in-plane but normal to 
the chain) axis, and present explicit expressions for the corre¬ 
sponding effective Hamiltonian. We note that due to the pres¬ 
ence of SOC, the effective Hamiltonian is anisotropic which 
can be readily seen already at the single-impurity level [62]. 

In the a) and b) cases, by transforming the effective tight- 
binding Hamiltonian 1-Ls{hj) to momentum space, we find 
that the corresponding Bogoliubov-de Gennes (BdG) Hamil¬ 
tonian reads 


To order and the effective hopping energy is 

given by 

/i,(fc)~e, + i[/o.+(fc) + /o.-(fc)] (18) 

hP + h^pk\ 

hx{k) €x + - [/o,+ (fc) + Jo,-(A:)] (19) 

hP+h^pk^ 

where the functions /o,±(^) and l 2 ,±{k) are obtained by tak¬ 
ing the E —> 0 limit in Eqs.(A17). The on-site energy is 

— Jo Q(^Ji(2 — Jo + Ji) 

“ Jo Jo(l + Jl)2 ’ 

( 20 ) 

1 - Jo [2(1 - Jo) + (1 - J2)] 

"" “ Jo Jo(i - j?y 

The effective p-wave pairing take the form 

A,(fc) Ax{k) [K,,+ {k) - K,,-{k)] 

+ ( 21 ) 

I Jl 

( 22 ) 

The functions In,±{k) = In,±{k,E = 0) and Kn^±{k) = 
Kn,±{k,E = 0) are defined in the Appendix A, see 

Eqs.(A17-A20). The functions h{k) and A{k) have the fol¬ 
lowing properties h{k) = h{—k) and A(—/c) = —A{k), and, 
thus, the gap is generically vanishing at k = 0,7r/a. More¬ 
over, the effective pairing A{k) is vanishing for a ^ 0 analo¬ 
gous to the semiconductor nanowire proposal [17, 18] in that 
the SOC controls the excitation gap. 

To have a better understanding of the Hamiltonian struc¬ 
ture in Eq. (17), it is instructive to perform a perturbative ex¬ 
pansion of h{k) and A{k), for example, around k = 0, see 
Eqs. (18) and (21) where the expressions for and 
are given in the Appendix E. In the limit a/^o ^ 0, where 
^0 = t’F/A is the superconducting coherence length, these 
functions have singular points for some values of kpct which 
is a consequence of the long-range nature of the hopping ma¬ 
trix element in the effective Hamiltonian. The presence of a 
finite coherence length however, regularizes the singulari¬ 
ties. Nevertheless, such a strong dependence on kpa leads to 
significant variations of the effective mass and Eermi veloc¬ 
ity. Another important feature that we find is that the effective 
Hamiltonian is anisotropic due to the SOC (cf. Eq.(20)) which 
might be helpful to drive the topological transition by chang¬ 
ing the direction of the magnetization of the impurities form¬ 
ing the chain. We note that this effect is absent for Ji = 0, 
in which case we recover the results of Ref. [43]. Thus, the 
dependence of the effective chemical potential on the angle 0, 
which is the only tuning parameter in the Hamiltonian (17), is 
a feature of the multichannel magnetic impurity model. 
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Finally, we discuss the case in which S || y. In this case the 
magnetic exchange energy term and the SOC term commute, 
and, as a result, the system is gapless. Thus, this case is not 
interesting from the point of view of Majorana physics. 


B. Topological Properties 



FIG. 2. (Color online) Topological phase diagram for the deep s 
band as a function of physical parameters. The dark and light colors 
represent topologically nontrivial and trivial phases, respectively, (a) 
Topological phase diagram for the magnetization in z direction, as a 
function of Jo and kpa for a = 0.3, Ji = 0.4, and = 2a. (b) 
Topological phase diagram for the magnetization in z direction as a 
function of a and kpa for Jq = 1.025, Ji = 0.4, and = 2a. 
(c) The phase boundary for the magnetization in z direction (blue 
dashed line) and for the magnetization in x direction (red solid line) 
indicates that by changing the magnetization one can drive the topo¬ 
logical phase transition, (d) Calculated quasiparticle excitation gap 
for the parameter regime in the phase diagram (a) with the phase 
boundary indicated by white line. 


Having derived the effective Hamiltonian (17), we can now 
calculate the topological phase diagram. The Hamiltonian 
(17) for a generic direction of magnetization is in the sym¬ 
metry class D [67-69], and, thus, is characterized by the Z 2 
topological invariant, the so-called Majorana number A4 [70]: 

j\4 = sgn [h{{))h(Ti/a )]. (23) 

The system is in the topological superconducting phase when 
M = —1, whereas A4 = +1 indicates a non-topological 
phase. We obtain the topological phase diagram by calculat¬ 
ing M. 

Figure 2 (a) shows the topological phase diagram in the 
{kFd^Jo) plane for the deep s-band limit for the case in which 


the magnetic moments of the impurities forming the chain are 
aligned along the z direction and a = 0.3. The range of values 
of kpCi has been chosen so that the inequality X/^/kpa 1, 
on which the treatment of the previous section relies, is well 
satisfied. From Fig. 2 (a) we see that, for a = 0.3 there is a 
large fraction of the (kpa^Jo) in which the chain is expected 
to be in a topological phase characterized by odd number of 
Majoranas at its ends. It is interesting to ask how the topologi¬ 
cal phase diagram is affected by the strength of the SOC. This 
question is addressed by the results presented in Fig. 2 (b) 
that shows the dependence of M on a and kpd for a fixed 
value of Jo . Again we notice that there is a large fraction of 
the (a, kpo) space in which the chain is expected to be in a 
topological phase. 

Experimentally it can be challenging to vary in a controlled 
way parameters such as a, Jq, and kpd and therefore to ver¬ 
ify the theoretical predictions shown in Fig. 2 (a), (b). How¬ 
ever, our multichannel treatment, contrary to the single chan¬ 
nel treatment [39, 43, 48], reveals that, the topological charac¬ 
ter of the chain state also depends on the direction of the mag¬ 
netization. This is shown in Fig. 2 (c) in which we can observe 
that the boundaries of the topological phase in the Jq) 
plane are different depending on the direction, z or x, of the 
magnetic moment of the impurities forming the chain. This 
result can be easily understood by considering that the on¬ 
site energy, Eqs.(20), depends on the direction of the impurity 
magnetization. The dependence of the topological index on 
the direction of the chain magnetization is very important be¬ 
cause it allows, in principle, to tune the chain in and out of a 
topologically nontrivial phase by varying a quantity that can 
be tuned and controlled experimentally by using an external 
magnetic field. The appearance, disappearance, of a zero bias 
peak as a function of the direction of the magnetization of the 
chain according to theoretical predictions like the ones pre¬ 
sented in Eig. 2 (c) would provide compelling evidence of the 
Majorana character of the observed zero energy states. 

In addition to the topological index (Majorana number), we 
have also calculated quasiparticle excitation gap as a function 
of Jo, a, and kpd, see Eig. 2(d) and Eig. 3(a)-(c). One can 
notice that the closing of the gap is consistent with the phase 
diagram shown in Eig. 2(a)-(b). Additionally, one has infor¬ 
mation regarding the magnitude of the gap deep in the topo¬ 
logical phase which is crucial for understanding the stability 
of the topological phase. Eig. 3(c) shows there also exist gap 
closing points inside the topological phase, not determined by 
calculating the topological index. To identify the location of 
these gap closing points, we plot the hopping energy spectrum 
the effective pairing energy spectrum lmA{k), and the 
energy spectrum of the particle band as a function of k and a 
in Eig. 3(d)-(f), respectively. Clearly, for certain values of a 
both hopping and pairing energies vanish at certain momen¬ 
tum points between 0 and tt resulting in the gap closing. In the 
limit Ji ^ 0 it can be shown analytically that these accidental 
zeros occur when the kp^F are commensurate with each other, 
i.e. when a = n7r/{kFa), with n G N. We will discuss the 
origin of these accidental gap closing points in more details in 
the next section. 
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FIG. 3. (Color online) Quasiparticle excitation gap Eg along different line cuts on the phase diagrams: Panel (a) along the line-cut A in 
Fig. 2(a) as a function of Jo at = 37.27r. Panel (b) along the line-cut B in Fig. 2(a) as a function of at Jo = 1.015. Panel (c) along 
the line-cut C in Fig. 2(b) as a function of a at /cfu = SG.Stt. Here the shaded area in (a)-(c) indicates the topologically nontrivial phase as 
shown in Fig. 2(a)-(b). To see the location of the gap closing points in Panel (c), we plot the hopping energy spectrum /i(/c) in panel (d), the 
effective pairing energy spectrum lmA{k) in panel (e), and the whole energy spectrum of the particle band in panel (f), as a function of k and 
a. 


IV. DEEP P BAND LIMIT 
A. Derivation of the Effective Hamiltonian 

We now discuss the deep p-band limit assuming that the 
energy of the / = ±1 states is lower than that of / = 0, i.e., 
Jo ^ Ji‘ In the limit Ji ^ 1 and a‘^Jo 1, the on-site en¬ 
ergy e « A(1 — is close to the midgap [62]. Once 

again, we assume that the / = ±1 states are well separated 
from the / = 0 state, and, therefore, I = 0 states can be inte¬ 
grated out. As a result, we obtain a tight-binding description 
for the p-bands with the s channel taken into account pertur- 
batively by allowing for the transitions through intermediate 
virtual I = 0 states. Note that there is a significant differ¬ 
ence with respect to the calculation in Sec. Ill since there are 
now two low-energy p-bands. Following the same procedure 
as in the previous section, we obtain the 8 dimensional matrix 
equation for the deep p-band limit: 

=0, (24) 

3 

where -i? i) is the 8 dimensional spinor for 

the p-channel states. The derivation of the matrix is 

presented in the Appendix F. Assuming that kpct ^ 1 and 
e —> 0, the p-bands have narrow bandwidth with the cen¬ 
ter of the bands being close to the midgap. One can then 
linearize Eqs.(24) with respect to E: 'M.p{E) ^ — 

• E, and neglect the energy dependence of the inter-site 
matrix lim£;^o by dropping the terms 

O and O with a 1 (i.e. we will keep 


henceforth the terms only up to 0{a‘^) and 0{a/^/k^).). 
After some algebra (see Appendix F for details), Eq. (24) can 
be written as 

= (25) 

3 

where (0). We then project onto 

the local basis of YSR states: ( 0i,+, 02,+, 02,-)^ 

where 0 i( 2 ),± are the eigen spinors of the single-impurity 
bound states with energy ±ei( 2 ) in channel 1 (/ = —1) and 
channel 2 (/ = 1). The local basis can be found by solv¬ 
ing the single-site equation = 0 as a special case 

of Eq. (24), where the bound state energies are determined 
from Det [M.p{E)] = 0, see Appendix G. After the projec¬ 
tion onto the local basis we obtain the effective Hamiltonian 
'Hp{k) describing the two coupled bands of the YSR chain in 
the deep-p band limit. 

Simple analytical expressions for 'Hp{k) can be obtained 
when the impurity spins are polarized normal to the plane, i.e. 
along z direction. In the case that the impurities are polarized 
in z direction, transforming the effective tight-binding Hamil¬ 
tonian 1-Lp{i^j) to momentum space, we obtain a two-band 
BdG Hamiltonian, 

/ hii{k) hi2{k) Aii(fc) Ai2(fc) \ 

h2l{k) h22{k) A2l{k) A22{k) 

A AJi(fc) A^i(fc) -hu{k) -h2i{k) ■ ^ 

\AUk) Al^{k) -/li2(fc) -h22{k) ) 

The coefficients here satisfy the following properties: 
hij{k) = hij{—k), Aij{—k) = —Aij{k) and, therefore. 
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FIG. 4. (Color online) The dependence of the normal-state energy 
spectrum (i.e. Aij = 0) on momentum k for (a) aJo = 0, Ji = 
1.0125, krci = ST.Stt, (fo = 2a; and (b) a = 0.3, Jo = 0.4, Ji = 
1.0125, kpci = 37.57r, = 2a. The normal-state spectrum consists 
of heavy-fermion and light-fermion bands which are hybridized by 
the SOC. The zoom-in figure of panel (b) near the Fermi level is 
shown in the inset. 


^ij{p) = 0 atj 9 = 0, 7 r/a. The effective hopping energies 
include both intra-channel hopping 


hii{k) ~ ei + - [Io,+ {k) + Io-{k )], 


(27) 


/i22(fc)2^e2 + -[Jo,+ (fc)+/o.-(fc)], 


(28) 


and inter-channel hopping 


huik) 


/l2l(fc)=i^[/2,+ (fc) + /2,-(fc)], 


(29) 


and inter-channel pairing 


Ai2(fc) = A2i(fc) ~ - [Ki, + (fc) - 


(34) 


ICxJ Q 


2(1 

A 


12 


Jo) 

k. 




The coefficients of the small k expansions are explained in the 
Appendix E. 

In order to understand the physics described by the Hamil¬ 
tonian (26), we first discuss the effect of SOC on the normal- 
state band structure (i.e. Aij = 0). The spectrum for the two 
bands reads 

hii{k) + h22{k) ± + {Sei2f 

(35) 

where dei 2 = ci — 62 . As shown in the Appendix A, to leading 
order in Xjyjkpa, Io,x{k) ^ l 2 ,x{k). Hence hi 2 is approxi¬ 
mately the same as hn and /i 22 - In the absence of SOC, 5ei2 
vanishes, and the band structure is characterized by a heavy- 
fermion band (1 — Ji)/Ji crossing with a dispersive 

light-fermion band ^ {1 — Ji)/Ji + [/o,+ (/c) + /o,-(/c)] 
with the bandwidth doubled compared to s-band, as shown in 
Fig. 4(a). The on-site orbital structure of these two bands are 
symmetric(light) and anti-symmetric (heavy) combinations of 
I = ±1 states. The physical origin of these orbital structures 
reflects the degeneracy due to the isotropic magnetic potential 
and the asymptotically equal hopping amplitudes. 

In the presence of SOC and a flnite s channel coupling (i.e. 
Jo 7 ^ 0 and a 7 ^ X)) 5ei 2 becomes flnite, which induces a hy¬ 
bridization gap between the two bands leading to an avoided 
level crossing, as shown in Fig. 4(b). The induced hybridiza¬ 
tion gap gives rise to an interesting feature in the topological 
phase diagram as discussed below. 


Em _ 1 

A 2 


with the on-site energies 


B. Topological Properties 


^ 1 — Ji Jq{2 — Ji + Jo) 
“ Jl Jl(l +Jo)2 
1-Jl 


(30) 

(31) 


The effective p-wave pairing also contains both intra-channel 
pairing 


All(fc)2.-[Kl,+ (fc)-i^l,_(fc)] 


(32) 


icxJo 

X Jo 

^11 


[2Ai,+ (A:)+2Ai,_(A:)], 




mk) ^ ^[Km) - Em)] 


(33) 




The topological phase diagram in the deep p-band limit 
(26) involves two bands which are hybridized by the SOC. 
As a result, it exhibits a more intricate dependence on the pa¬ 
rameters compared with the deep 5 -band limit. In order to 
calculate the Z 2 topological invariant A4, we use the method 
developed for the multiband semiconductor nanowire system 
[71]: M = sgn [Pf 5 ( 0 )Pf 5 ( 7 r/a)] = ± 1 , where the an¬ 
tisymmetric matrix B{p) = l-Lp{p)Tx is calculated at the 
particle-hole invariant points: p = 0, n/a. For the two- 
band system, the corresponding expression for the Pfafflan is 
PfB(p) = hi2{p)h2i{p) - hii{p)h22{p)- 
The topological phase diagram for the effective Hamilto¬ 
nian (26) as a function of Ji and kpci is shown in Fig. 5 (a). 
We immediately notice that in the deep p-band limit there is 
a large fraction of the (Ji, plane in which the chain is 

in a topological phase. This is somewhat surprising given that 
in the deep p-band we have two bands close in energy and 
one could expect that the hybridization of the two bands could 
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FIG. 5. (Color online) (a) Topological phase diagram in the 
{kpa, Ji ) plane in the deep p-band limit, for the case in which S || z, 
Jo = 0.4, a = 0.3, and = 2a. (b) Enlargement of the topolog¬ 
ical phase diagram shown in (a) around the region surrounded by 
the dashed line rectangle, (c) Calculated quasiparticle excitation gap 
for the parameter regime in the phase diagram (b) with the phase 
boundary indicated by white line, (d) The quasiparticle excitation 
gap on the line-cut A in panel (b) and (c) near the re-entrance region 
at kpa = 36.47r. (e) The quasiparticle excitation gap on the line-cut 
B in panel (b) and (c) near the re-entrance region at kpa = 37.67r. 
Here the shaded area in (d) and (e) indicates the topologically non¬ 
trivial phase as shown in panel (b). The integer numbers -1, 0, 1 
shown in (d) and (e) are the winding numbers calculated for each 
gapped phase. 


lead for most values of Ji and kpa to a situation in which the 
Fermi energy intersects the bands an even number of times. 
The reason why this is not the case is evident from the plots 
of the two bands for a typical situation shown in Fig. 4: one 
can see that the result of the hybridization of the bands leads 
to the formation of the light- and heavy-fermion bands. As a 
consequence, for most values of Ji and the Fermi energy 
crosses only once the light-fermion band. On the other hand, 
we expect that when Ji and kpa are such that the Fermi en¬ 
ergy is very close to the heavy band the chain should be in 
a topologically trivial phase. Given the flatness of the heavy 
band we expect that this will happen only for a very small 
range of values of Ji. Indeed, a close inspection of the topo¬ 
logical phase diagram of Fig. 5 (a) shows that for Ji very 
close to 1 there is a narrow region, highlighted in the figure 
by a rectangular box, in which the chain is in a topologically 
trivial phase. 

To understand the re-entrance region from the topological 
phase to the trivial phase when Ji ~ 1, in Fig. 5 (b) we show 


a zoom-in of the topological phase diagram in the region sur¬ 
rounded by the rectangular box in Fig. 5 (a). In the zoom-in 
plot the presence of a topologically trivial phase is clearly vis¬ 
ible. In this region, due the fact that Ji ^ 1 the effective 
chemical potential is close to zero so that the heavy-fermion 
band becomes important: in this situation the chemical poten¬ 
tial lies inside the hybridization gap, as shown in Fig. 4 (b), 
which results in the re-entrance to the trivial phase. Outside 
this narrow range of values of Ji the heavy-fermion band is 
either completely empty or filled and the chemical potential 
crosses the light-fermion band an odd number of times, thus, 
there is an odd number of Majorana zero modes per end orig¬ 
inating from the light-fermion band. Due to this interplay, the 
two YSR bands do not annihilate each other (except for the 
small region close to Ji ^ 1) which is a peculiar feature of 
this model. 

In order to understand the stability of the topological phase 
and to corroborate the results discussed above, we also com¬ 
pute the quasiparticle excitation gap Eg, see Fig. 5(c). Fig¬ 
ure 5(d)-(e) plot the value of Eg along the two line-cuts on the 
phase diagram near the re-entrance region. One can see that 
the quasiparticle gap closing is consistent with the topologi¬ 
cal phase diagram in Fig. 5(b). We once again find that SOC 
controls the magnitude of the quasiparticle gap and, as such, 
is crucial for the stability of the topological phase. 

In addition to the quasiparticle gap closing at the topolog¬ 
ical phase boundary, one can notice that there are also points 
where the gap vanishes in the topological phase, see Fig. 5(d) 
at Ji ^ 1.005 and Fig. 5(e) at Ji ^ 1.017. We now inves¬ 
tigate in detail the two regions across these points using an 
additional symmetry of our effective Hamiltonian (26). In 
addition to the particle-hole symmetry V = TxE where JC 
refers to complex conjugation, our effective spinless Hamilto¬ 
nian also has a pseudo-time reversal symmetry T = 1C. Using 
these two symmetries, one can construct another symmetry 
- chiral symmetry S = TV = Tx which anticommutes with 
the Hamiltonian (26). Thus, the Hamiltonian (26) belongs 
to the BDI symmetry class [67-69] which is characterized 
by the integer invariant W and supports multiple spatially- 
overlapping Majorana zero modes [72]. In order to calcu¬ 
late the topological index W, it is convenient to transform the 
Hamiltonian (26) into a chirality basis using a unitary trans¬ 
formation U = which converts the Hamiltonian to the 

off-diagonal form: 

U%;{k)U^ = ^0^^) ■ (36) 

Then, the winding number (i.e. the number of Majorana 
modes per each end) can be calculated by introducing a com¬ 
plex variable z{k) = det[A{k)]/\ det[A(/c)] |, and calculating 
the integral 


>v = -i 

TT Uo Z{k) ’ 


(37) 


Using this analysis we find that, for example, the phases at 
Ji = 1.001 and Ji = 1.007 in Fig. 5(d) have different wind¬ 
ing numbers >V(Ji = 1.001) = 1 and>V(Ji = 1.007) = —1. 
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Thus, gap closing between these two regions corresponds to 
the transition between W = ±1. The same argument holds 
for Ji = 1.015. Thus, accidental gap closing points inside of 
the topological or non-topological phases are not really acci¬ 
dental but represent the change of the winding number by an 
even integer. 

The analysis above relies on the chiral symmetry. However, 
in realistic systems the chiral symmetry can be easily broken 
by allowing, for example, for a generic direction of magnetic 
chain polarization (i.e. along ^-axis). The precise magnitude 
for the Majorana splitting energy, which is important for tun¬ 
neling transport measurements, depends on the details of the 
chiral-symmetry-breaking perturbations, and we refer a reader 
to Refs. [45, 48, 73-75] for more details. As a consequence, 
the topological phases identified by the parity of the topolog¬ 
ical index are expected to be much more robust and this is the 
reason that our analysis has been focused mostly on character¬ 
izing the dependence of such index on the parameters of the 
system. 


V. CONCLUSIONS 

We have studied the topological properties of a chain of 
magnetic impurities placed on the surface of an s-wave su¬ 
perconductor with Rashba spin-orbit coupling taking into ac¬ 
count the presence of multiple scattering channels, in the limit 
in which the states induced by isolated impurities are well 
described as Yu-Shiba-Rusinov states and the distance a be¬ 
tween the impurities forming the chain is such that kpa ^ 1. 
The inclusion of multiple angular momentum scattering chan¬ 
nels, /, implies that for the chain of YSR states we have multi¬ 
ple bands. We have shown that the multiband character of the 
bands strongly affects the topological properties of the chain 
for the case when the lowest energy bands are the ones arising 
from the hybridization of YSR states with |/| > 0. Consider¬ 
ing the lowest I = 0 and | /1 = 1 channels we have obtained the 
topological phase diagram in the deep 5-band and deep p-band 
limits. Our results show, somehow unexpectedly, that even in 
the deep p-band limit there is a large region of parameters 
space in which the chain is in a topological phase. Moreover, 
we find that even though the deep p-band case involves two 
bands, and it is not obvious a priori that these two bands do 
not “annihilate” each other, this limit seems to be more favor¬ 
able for the observation of Majorana zero modes. This can be 
seen from the comparison. Fig. 6, of the topological phase di¬ 
agram obtained in the deep 5-band and deep p-band limit: we 
see that, for the same range of values of kpa and relevant cou¬ 
pling constants (Jo in the deep 5-band, Ji in the deep p-band 
limit) the phase space where the chain is in a topologically 
non-trivial state is larger in the deep p-band limit than in the 
deep 5-band limit. We have also characterized the stability of 
the topological states by computing the quasiparticle excita¬ 
tion gap. 

I 



(b)u 


Ji 



rp 


kfa/n 


36.0 36.5 37.0 37.5 

kfa/n 


FIG. 6. (Color online) Comparison of the topological phase diagram 
in (a) the s band system for o = 0.3, Ji = 0.4, and = 2a; and in 
(b) the p band system for a = 0.3, Jo = 0.4, and = 2a. 


Our results have important implications for the ongoing ex¬ 
perimental search for Majorana zero modes in this system. 
The interatomic spacing in Ref. [52] is of the order of the 
Fermi wave length (kpa ^ l),m which case direct tunneling 
between iron atoms needs to be included, whereas our calcu¬ 
lation assumes kpa ^ 1. However, at the qualitative level 
one can already draw a number of conclusions. As kpa de¬ 
creases the bandwidth of subgap YSR bands becomes larger 
and, thus, the mixing between them becomes even more im¬ 
portant. Therefore, it is conceivable that the subgap features 
reported recently by Yazdani et al. [76] might correspond to 
/ > 0 YSR bands. In order to understand these additional sub¬ 
gap features within the framework of our model, one would 
need to systematically measure the spectrum of the YSR states 
of single- and two-atom structures [62] to identify the angu¬ 
lar momentum channels that are the most relevant for the one¬ 
dimensional chains. Another important feature that we predict 
is that SOC allows one to tune the effective chemical potential 
for the YSR bands and drive the topological phase transition 
using an external magnetic field. As such, large SOC might 
be helpful for the manipulation of the Majorana zero modes 
in this system. At the theoretical level, it would be interesting 
to establish the correspondence between our results involv¬ 
ing classical magnetic impurities with the more microscopic 
multi-orbital Anderson model for magnetic impurities. 
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Appendix A: Analytic Expressions of Integrals a and a 


In this Appendix, we provide analytic expressions of the integral functions defined in Eq. (11) and Eq. (12): 

IoA^-,E) = -0^Re [Jo {{kF,x+iCx^)\x\)+iHo + ^Ca') kl)] 

Ko,x{x; E) = 7 Alm [Jq {{kF,x + Ka”')!^!) + ((^f.a + iCx^)\x\)] 


h,\{x;E) = -sgn[a;] 


iAjx 


a/A 2 - j ;2 


Re [Ji ((fcF.A +<A')kl) - (fe +<A')kl)] 


Ki^x(x\E) = sgn[a:;]i7Alm [Ji ((fc^.A + <a ^)|2^I) “ ((^f,a + Kx ^)kl)] 

J 2 {{kF,x + Kx^)\x\)+iH _2 {{kF,x+iCx^)\x\) + 


l2,x{x;E) = -. .^'^^ _^ Re 


VA 2 - ^2 


2i 


K 2 ,x{x;E) = - 7 a Im 


J 2 {{kF,\+iCx ^)kl) +'iH-2 {{kF,\+iCx ^)kl) 


TT {kF,X + iCx kl 

2i 


h^x{x\E) = sgn[a;] _ Re 


a/A 2 - £’2 


n {kF,x + iCx 

J 3 {{kF,x + iCx^)\x\) -iH-3 ((^f.a + <a^)| 2 ^I) + 


6i 


K3,x{x\E) = -sgn[a:;]i 7 Alm 


J 3 ((^F.A + Kx ^)kl) - {{kF,X + Kx ^)|2:|) + 


F [(A:f,a+<a^) kl] . 

6 i 


F [(A:f,a+<A^) 1 ^ 1 ] . 


(Al) 

(A2) 

(A3) 

(A4) 

(A5) 

(A6) 

(A7) 

(A8) 


Here Jn{z) and Hn{z) are Bessel and Struve functions of order n, respectively; Ca ^ , and 7 a = 1 + A ^^“^2 • Note 

that the expressions for Ki^x{x', E) given above are valid for a; ^ 0, and the integral Ki^xiK E) = 0 for x = 0. Assuming 
fcF|a:| 1 and Ca ^ <C kF,x, we can use the asymptotic forms of the Bessel and Struve functions [77], In the limit 

kFX » 1, one can find approximate expressions up to the order l/{kFx)‘^: 

-A7a 


Io,x{x; E) = 


\/A2 - j;2 Y T^kF,x\x 


^ -e 


cos{kF,x\x\ - If) + ^^-*^sin(A:F,A|a;| - ^f) 


Ko,x{^’E) =7ai 


Ii^x{x;E) = -sgn[a:] 


FfcF.Akl 
iA-yx 


e 


sin(fcF,Akl - ^f) - cos(fcF,Akl - ^f) 


27a 


nkF,x\x\ 


a/A2 - E^]I TTkF,x\x 

Ei,x{^'’ E) = sgn[a;]t7AW ^ 


2 


cos(fcF,A|a:| - ^f) - sin(fcF,A|a:| - ^f) 


(A9) 

(AlO) 

(All) 


irk 


F,xm 


33 3 

sin(fcF,A|a:| - -f) + cos(fcF,A|a;| - -f) 


+ sgn[a;] 


*27a 


F(A:F,A|a:|) 

(A12) 


l 2 x{x;E) = —=^2^. 

^ VA2 - E^ 




2 

e 




'Cx l^l 


sin( 


1 15 1 

cos(A:F,A|a:| - -f) ~ sin(fcF,Ak| - ^f) 

27a 


F^F.Ala;, 


■f) 


FfcF.Akl 




c»s{k,M - 


-.A*; E) = ■I-I [.i„,fc„|*| ^ 


K- 


TTkF,x\x\ 

cos(fcF,; 


• tt- I I 3 . 

^ - j,) 


ki - 


r 1 ^ 

-sgn[x]^— 

TT [ki 
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The corresponding Fourier transforms of the above asymptotic forms to the leading order of 77 ^^ are given by 


IoAk;E) = l2Ak'.E) = 


e“*3^Lii +e*5^Lii (e“*' 


-A 7 A / 1 

a/A 2 - £’2 y 27rA:F,Aa 


—ao-— C\ 


ika^ 


—ikF,xa—Cx CL— 


^a—ika^ 


(A17) 


-S') = K2,\{k E) = -ijx] 


1 


^irkp xa 


e '4 


J^Lii (e 




z/cf^ag-Ca ) — e^^^Lii ( 


+ e“*4^Lii 


^ikF,xa-Cx^ci-ika^ ^-ikF,xCL-Cx^CL-ika^ 


(A18) 


/i,a(A:;F;) = / 3 ,a(A:;F;) = 


iA7A 


e“*I^Lii + e*?^Lii (^-'^^F,xa-Cx^a+ika^ 

(A19) 


VA2 - E‘^]J 27rkF,xa 

— ^^^F,xa-Cx^a-ika^ — e^i^lAi ^-'^^F,xa-Cx^a-ika^ 


K^^x{hE)=K^^x{hE)=-fx^ 


1 


27r k F xCi 


e“*f^Lii ^^^^F,xa-Cx^CLFika'^ — e^^^Lii ^^-'^f^F,xCi-Cx^CLFika'^ 


-e“^4^Li 


ii + e^^^Lii ^^-^^F,xa-Cx^ci-ika'^ 


where Lis( 2 ;) is the polylogarithm function 


Li«(^) = E;;7- 


n=l 


(A20) 


Appendix B: Calculation of the Greens functions {E) 

The Green’s function of a superconductor with Rashba spin-orbit coupling can be expanded in the angular momentum channels 


G(k; E) = Gi{k; The local Green’s function (i.e. i = j) reads 

/ hfjh _ 

—Gi.i,(k-,E)^Gi.i,{E). (Bl) 

Using Eqs. (4)-(6), one finds that there are three non-zero local Green’s function corresponding to angular momenta / = — 1, 0, 1: 

G-iiE) = ^ 2^2 + Act+t^) , (B2) 

Go{E) = ^EaoTo + Actot,) , (B3) 

Gi(E) = -^||^=^(E<T_ro + A<T_r,). (B4) 
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The non-local Green’s function (i ^ j) is defined as 




A=ili 


E 


— Cr+To + CF-^Tx ) /|_i+^|,A + 


+ 


E 


-^O-qTq + (JQ^x ) ^1 m| ,A + (cTor^) X| m| ,A 


-h (iA) 




— Cr_ro H- (J-Tx 1 /|l+m|,A + ^|l+m|,A 


(B5) 


where m = / — /' = 0, ±1, ±2, and the functions In,x{xij; E) and Kn^x{xij]E) are given in Appendix A. 

Appendix C: Derivation of effective Hamiltonian in deep s-band limit 

In this Appendix we provide the details of the derivation of effective Hamiltonian in deep 5-band limit. Using the method 
outlined in Sec. II, we find that the effective eigenvalue equation for the deep 5-band limit reads 

= (Cl) 

hi 

After substituting Eq. (14) back into the above equation, one finds 

0 =(M“o - q + j) 

+ (C2) 

The assumption kpct ^ I allows one to neglect the terms Using the tight-binding approximation, one finally 

arrives at 


0 = 




'•Pifi 


-E 

-E 

3^^ ^ 


^^•0 + 0 {{Mh^) 


P’i 


3,0 




(C3) 


Appendix D: Local basis in the deep s-band limit 

In this Appendix, we provide details for the projection procedure used to the derive an effective Hamiltonian for the deep 
5-band limit. Assuming that the magnetic-atom polarization is along z axis, the unnormalized local basis for magnetic impurity 
reads: 

^ ( 1, 0, 1, 0 and ^ ( 0, 1, 0, —1 . 

When the polarization is along x-axis, the unnormalized local basis for magnetic impurity is given by 

</.+ -(!, 1, 1, l)^and^_~(-l, 1, 1, -1)^. 
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For the sake of completeness, we also derive effective Hamiltonian when the magnetization is along ^-axis. In this case, the 
unnormalized local basis for magnetic impurity becomes 

^ ( 1 , i, 1 , i and (^_ ^ ( 1 , ——1 . 

After the projection onto the local basis, the effective Hamiltonian is given by 


um 


hy{k) + dy{k) 0 

0 —hy{k)+dy{k) 


(Dl) 


with the functions hy{k) and dy{k) being 


hy{k) — ey + - [Io^+{k) + /o,_(fc)], 




J\OL 


(1 - Jl) 




(D2) 

(D3) 


Here the on-site energy Cy = Cx. One can see that quasiparticle spectrum in this case is indeed gapless. 


Appendix E: Effective Hamiltonian in the long wavelength limit k ^ 0 


It is instructive to expand the functions I (n, k) and K (n, k) appearing in our effective Hamiltonian close to /c = 0 in order to 
understand the spectrum qualitatively. After some algebra, one finds 


Io,x{k; E = 0)= l2,x{k; E = 0) = 
h,x{k; E = 0)= h,x{k\E = 0) = - 27 ^^ 
E = Q)= K3^xik;E = 0) = 2jx^ 

where 


1 


I 1 

irkp xa 


Ao{kF,xa + ^Ca ~ 27a 


1 


Trkpx^ 


Bi{kF,xCi + iCx ): 


1 


irkF^x^^ 


Ci{kF,xa^iCx ci)k^O{k^), 


Ao{z) = Re [lIi (e^^)j + Im \Ui (e^^) 

A2{z) = Im Li_3 (e*^) 

Bi{z) = Re Li_i (e*^) — Im Li_i (e*^) 
Ci{z) = Re [Li_i (e^^)j + Im [Li_i (e^^) 


A2{kF,xa + Kx + C’(fctfil) 

(E2) 

(E3) 

(E4) 

(E5) 

(E 6 ) 

(E7) 


The dependence of the functions ^ 0 ( 2 ;), ^ 2 ( 2 ;), Bi{z) and Ci{z) on the external parameters is shown in Fig. 7. One can 
notice that when kF^x^ = 27rn with n being an integer, these functions have singularities which follows from the definition 
of polylogarithm function. These singularities are cutoff by the finite coherence length. In realistic systems, however, the 
superconducting coherence length is much larger than the interatomic spacing, and, thus, the parameters such as effective mass 
and Fermi velocity are strongly dependent on see Fig. 7. 
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FIG. 7. (Color online) The dependence of the functions Ao(kF,xCi-^iCx A2{kF,xCi-^iCx Bi(kF,xCi-^iCx Ci{kF,xCi-^iCx 

on kF,xa- Here we used Ca = 10a. 
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Finally, the expansion of the coefficients in the deep 5-band Hamiltonian at /c 0 becomes 


C='.- (‘+^'1 


*'‘’=E‘TA\/rf^ +.<,-■«). 

4 '’=-ET>\/j^ (^ + +«.-■»)■ 


The expansion coefficients in the deep p-band Hamiltonian are 


=ei ~r — Ao{kF,\a + iC^^a), 

^ y TTkF,\a 

>^22 = ^2 - y^TA-iZ-TT—^o(A:F,Aa + <A^“)’ 

A V 

= = - '^'y>'\ -r- — Ao{kF,\a + iC^'^a), 

^ y TTfcF.Aa 

= ^ 2 ^ = ^ 12 ^ = — A2{kF,\a + iCx^a), 

^ y TTKfxCL 


^n’ = E y - S) + •<»“»). 

^ 2 ^ ~i~ — Ci{kF,\a + iCx^a), 

^ y TrfcF.Afl 

= - CiXkF,xa + iC^^a). 


TTkF,\a V 1 + Jo 


(E8) 

(E9) 

(ElO) 

(Ell) 

(E12) 

(E13) 

(E14) 

(E15) 

(E16) 

(E17) 

(E18) 

(E19) 

(E20) 

(E21) 

(E22) 


Appendix F: Derivation of effective Hamiltonian in deep p-band limit 

In this Appendix, we provide the details of the derivation of effective Hamiltonian in deep p-band limit. Using the method 
outlined in Sec. II, we find the effective eigenvalue equation for the deep p-band limit. The corresponding equations for p-wave 
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bands are given by 


3,1 

3,1 

In order to integrate out s-channel, we have to solve for 1 p^^Q finding that 


(FI) 

(F2) 


3^i,l 


(F3) 


Substituting Eq. (F3) into Eq. (El, E2) and following the same procedure as in Appendix C, we eventually obtain two coupled 
equations for the p-wave bands 


-E 


-E 

37^^ 


0 = - + (M«i - M«o(M«o)-'M«i))A,,i 

mEi - mE(ME)-'mEi - ME(ME)-'mEi + ME(ME)-'mE(ME)-'mEi) 


'4’jA (F4) 


-E 


-E 


Ml, - + M1*o(M»o)-'M*^o(M“o)-iM“i 

After some manipulations, one can write eigenvalue equations in the compact form, see Eq. (24). 


V’i.i 


(F5) 


Appendix G: Local basis in the p-band system 

In this Appendix, we discuss the projection procedure in the p-band limit in the case of magnetic impurity spins being aligned 
along z axis. In this case, the unnormalized local spinors are given by 

</>!,+-( 1 , 0 , 1 , 0 , 0 , 0 , 0 , 0 )^, </. 2 ,+ ~( 0 , 0 , 0 , 0 , 1 , 0 , 1 , 0 )^; 


01,_~(O, 0, 0, 0, 0, 1, 0, -1)^, </.2,--(0, 1, 0, -1, 0, 0, 0, 0)^. 


Projecting onto the local basis ( 0i,+, 02,+, 01,-, 02,- ) ^, one arrives at the two-band Hamiltonian defined in Eq. (26)-(34). 
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